#include<iostream>
#include<vector>
#include<queue>
#include<malloc.h>
#include <math.h>
#include <time.h>

void R_K()

{
  int i,n,q;
  double u,k1,k2,T1,T2,z1,z2;
  double *v1,*v2,*v4,*v5,*u1,*u2,*u4,*u5,*y1,*y2,*y3,*y4;
  double e1,e2,emaxnorm,premaxnorm,Crate;
  double start,finish,duration;


  u1=(double *)malloc(sizeof(double)*10000000);
  u2=(double *)malloc(sizeof(double)*10000000);
  u4=(double *)malloc(sizeof(double)*10000000);
  u5=(double *)malloc(sizeof(double)*10000000);
  
  v1=(double *)malloc(sizeof(double)*10000000);
  v2=(double *)malloc(sizeof(double)*10000000);
  v4=(double *)malloc(sizeof(double)*10000000);
  v5=(double *)malloc(sizeof(double)*10000000);

  y1=(double *)malloc(sizeof(double)*10);
  y2=(double *)malloc(sizeof(double)*10);
  y3=(double *)malloc(sizeof(double)*10);
  y4=(double *)malloc(sizeof(double)*10);

  u1[0]=0.994;u2[0]=0;u4[0]=0;u5[0]=-2.0015851063790825224;u=1/81.45;
  T1=17.06521656015796;T2=19.14045706162071;


      std::cout<<"the classical RK method(situation1):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;
    
  //
 
  
  n=5e3;k1=T1/n;
  for (i=1;i<=n;i++){
    z1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);z1=pow(z1,1.5);
    z2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z2=pow(z2,1.5);

    y1[1]=u4[i-1];
    y1[2]=u5[i-1];
    y1[4]=2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/z1-((1-u)*(u1[i-1]+u))/z2;
    y1[5]=u2[i-1]-2*u4[i-1]-(u*u2[i-1])/z1-((1-u)*u2[i-1])/z2;
    
    v1[i-1]=u1[i-1]+k1*y1[1]/2;
    v2[i-1]=u2[i-1]+k1*y1[2]/2;
    v4[i-1]=u4[i-1]+k1*y1[4]/2;
    v5[i-1]=u5[i-1]+k1*y1[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y2[1]=v4[i-1];
    y2[2]=v5[i-1];
    y2[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y2[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y2[1]/2;
    v2[i-1]=u2[i-1]+k1*y2[2]/2;
    v4[i-1]=u4[i-1]+k1*y2[4]/2;
    v5[i-1]=u5[i-1]+k1*y2[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y3[1]=v4[i-1];
    y3[2]=v5[i-1];
    y3[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y3[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y3[1];
    v2[i-1]=u2[i-1]+k1*y3[2];
    v4[i-1]=u4[i-1]+k1*y3[4];
    v5[i-1]=u5[i-1]+k1*y3[5];

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y4[1]=v4[i-1];
    y4[2]=v5[i-1];
    y4[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y4[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    u1[i]=u1[i-1]+(k1/6.0)*(y1[1]+2*y2[1]+2*y3[1]+y4[1]);
    u2[i]=u2[i-1]+(k1/6.0)*(y1[2]+2*y2[2]+2*y3[2]+y4[2]);
    u4[i]=u4[i-1]+(k1/6.0)*(y1[4]+2*y2[4]+2*y3[4]+y4[4]);
    u5[i]=u5[i-1]+(k1/6.0)*(y1[5]+2*y2[5]+2*y3[5]+y4[5]);
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);

  //

  
  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;

  for (i=1;i<=n;i++){
    z1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);z1=pow(z1,1.5);
    z2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z2=pow(z2,1.5);

    y1[1]=u4[i-1];
    y1[2]=u5[i-1];
    y1[4]=2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/z1-((1-u)*(u1[i-1]+u))/z2;
    y1[5]=u2[i-1]-2*u4[i-1]-(u*u2[i-1])/z1-((1-u)*u2[i-1])/z2;
    
    v1[i-1]=u1[i-1]+k1*y1[1]/2;
    v2[i-1]=u2[i-1]+k1*y1[2]/2;
    v4[i-1]=u4[i-1]+k1*y1[4]/2;
    v5[i-1]=u5[i-1]+k1*y1[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y2[1]=v4[i-1];
    y2[2]=v5[i-1];
    y2[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y2[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y2[1]/2;
    v2[i-1]=u2[i-1]+k1*y2[2]/2;
    v4[i-1]=u4[i-1]+k1*y2[4]/2;
    v5[i-1]=u5[i-1]+k1*y2[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y3[1]=v4[i-1];
    y3[2]=v5[i-1];
    y3[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y3[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y3[1];
    v2[i-1]=u2[i-1]+k1*y3[2];
    v4[i-1]=u4[i-1]+k1*y3[4];
    v5[i-1]=u5[i-1]+k1*y3[5];

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y4[1]=v4[i-1];
    y4[2]=v5[i-1];
    y4[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y4[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    u1[i]=u1[i-1]+(k1/6.0)*(y1[1]+2*y2[1]+2*y3[1]+y4[1]);
    u2[i]=u2[i-1]+(k1/6.0)*(y1[2]+2*y2[2]+2*y3[2]+y4[2]);
    u4[i]=u4[i-1]+(k1/6.0)*(y1[4]+2*y2[4]+2*y3[4]+y4[4]);
    u5[i]=u5[i-1]+(k1/6.0)*(y1[5]+2*y2[5]+2*y3[5]+y4[5]);
  }

  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
  Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

      T1=T2,u1[0]=0.87978;u2[0]=0;u4[0]=0;u5[0]=-0.3797;
      std::cout<<"the classical RK method(situation2):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;

      n=100;k1=T1/n;
  for (i=1;i<=n;i++){
    z1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);z1=pow(z1,1.5);
    z2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z2=pow(z2,1.5);

    y1[1]=u4[i-1];
    y1[2]=u5[i-1];
    y1[4]=2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/z1-((1-u)*(u1[i-1]+u))/z2;
    y1[5]=u2[i-1]-2*u4[i-1]-(u*u2[i-1])/z1-((1-u)*u2[i-1])/z2;
    
    v1[i-1]=u1[i-1]+k1*y1[1]/2;
    v2[i-1]=u2[i-1]+k1*y1[2]/2;
    v4[i-1]=u4[i-1]+k1*y1[4]/2;
    v5[i-1]=u5[i-1]+k1*y1[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y2[1]=v4[i-1];
    y2[2]=v5[i-1];
    y2[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y2[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y2[1]/2;
    v2[i-1]=u2[i-1]+k1*y2[2]/2;
    v4[i-1]=u4[i-1]+k1*y2[4]/2;
    v5[i-1]=u5[i-1]+k1*y2[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y3[1]=v4[i-1];
    y3[2]=v5[i-1];
    y3[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y3[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y3[1];
    v2[i-1]=u2[i-1]+k1*y3[2];
    v4[i-1]=u4[i-1]+k1*y3[4];
    v5[i-1]=u5[i-1]+k1*y3[5];

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y4[1]=v4[i-1];
    y4[2]=v5[i-1];
    y4[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y4[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    u1[i]=u1[i-1]+(k1/6.0)*(y1[1]+2*y2[1]+2*y3[1]+y4[1]);
    u2[i]=u2[i-1]+(k1/6.0)*(y1[2]+2*y2[2]+2*y3[2]+y4[2]);
    u4[i]=u4[i-1]+(k1/6.0)*(y1[4]+2*y2[4]+2*y3[4]+y4[4]);
    u5[i]=u5[i-1]+(k1/6.0)*(y1[5]+2*y2[5]+2*y3[5]+y4[5]);
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);

  //

  
  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;

  for (i=1;i<=n;i++){
    z1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);z1=pow(z1,1.5);
    z2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z2=pow(z2,1.5);

    y1[1]=u4[i-1];
    y1[2]=u5[i-1];
    y1[4]=2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/z1-((1-u)*(u1[i-1]+u))/z2;
    y1[5]=u2[i-1]-2*u4[i-1]-(u*u2[i-1])/z1-((1-u)*u2[i-1])/z2;
    
    v1[i-1]=u1[i-1]+k1*y1[1]/2;
    v2[i-1]=u2[i-1]+k1*y1[2]/2;
    v4[i-1]=u4[i-1]+k1*y1[4]/2;
    v5[i-1]=u5[i-1]+k1*y1[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y2[1]=v4[i-1];
    y2[2]=v5[i-1];
    y2[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y2[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y2[1]/2;
    v2[i-1]=u2[i-1]+k1*y2[2]/2;
    v4[i-1]=u4[i-1]+k1*y2[4]/2;
    v5[i-1]=u5[i-1]+k1*y2[5]/2;
    

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y3[1]=v4[i-1];
    y3[2]=v5[i-1];
    y3[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y3[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    v1[i-1]=u1[i-1]+k1*y3[1];
    v2[i-1]=u2[i-1]+k1*y3[2];
    v4[i-1]=u4[i-1]+k1*y3[4];
    v5[i-1]=u5[i-1]+k1*y3[5];

    z1=v2[i-1]*v2[i-1]+(v1[i-1]+u-1)*(v1[i-1]+u-1);z1=pow(z1,1.5);
    z2=v2[i-1]*v2[i-1]+(v1[i-1]+u)*(v1[i-1]+u);z2=pow(z2,1.5);

    y4[1]=v4[i-1];
    y4[2]=v5[i-1];
    y4[4]=2*v5[i-1]+v1[i-1]-(u*(v1[i-1]+u-1))/z1-((1-u)*(v1[i-1]+u))/z2;
    y4[5]=v2[i-1]-2*v4[i-1]-(u*v2[i-1])/z1-((1-u)*v2[i-1])/z2;

    u1[i]=u1[i-1]+(k1/6.0)*(y1[1]+2*y2[1]+2*y3[1]+y4[1]);
    u2[i]=u2[i-1]+(k1/6.0)*(y1[2]+2*y2[2]+2*y3[2]+y4[2]);
    u4[i]=u4[i-1]+(k1/6.0)*(y1[4]+2*y2[4]+2*y3[4]+y4[4]);
    u5[i]=u5[i-1]+(k1/6.0)*(y1[5]+2*y2[5]+2*y3[5]+y4[5]);
  }

  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
  Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

}

    

    
    

    
    
    
